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~^With  on  exception,  all  finite  element  software  systems  have  element 
libraries  in  which  the  approximation  properties  of  elements  are  frozen.  The 
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J'Jrrs— syetaa,  called  COMET ^ploys  conforming  elements  based  on  complete  poly- 
nomials of  arbitrary  order.  The  elements  are  hierarchic^ i.e,  the  stiffness 
matrix  of  each  element  is  embedded  in  the  stiffness  matrices  of  all  higher  order 
elanents  of  the  same  kind.  The  user  controls  not  only  the  number  and 
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ttistribution  of  finite  elements  but  their  approximation  properties  as  well. 

Thus  convergence  can  be  achieved  on  fixed  mesh.  This  provides  for  very  effici- 
ent and  highly  accurate  approximation  and  a new  method  for  computing  stress 
intensity  factors  in  linear  eleastic  fracture  mechanics,  The  theoretical  devel 
raents  are  outlined,  numerical  examples  are  given  and  the  concept  of  an  advanced 
self-adaptive  finite  element  software  system  is  presented. 
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2.  CL’RRgNT  AND  AlfTICIPATED  RESULTS  Ai>mCIPATED  BY  THE  EtJD  OF  IHITIAL  GRANT 
PERIOD 

The  ultimate  objective  of  our  research  is  to  provide  algorithms,  and 
experimental  digital  computer  programs  which  implement  the  constraint  method 
for  three  dimensional  analysis,  emd  to  demonstrate  the  effectiveness  of  the 
method  in  three  dimensions.  Since  the  approaches  to  be  used  in  deriving 
these  algorithms  are  based  on  similar  approaches  for  two  dimensional  ema- 
lyses,  the  first  part  of  the  initial  grant  period  was  spent  in  completing  the 
formulation  of  algorithms  in  two  dimensions  so  as  to  have  a firm  foundation 
for  the  more  complex  three  dimensional  work.  Also  in  response  to  a suggestion 
by  Dr.  V.  3.  Venkayya  of  the  Air  Force  Flight  Dynamics  Laboratory  at  Wright- 
Patterson  AFB,  the  plate  bending  element  of  the  Constraint  Method  is  being 
developed  further.  The  results  obtained  and/or  anticipated  for  two  dimen- 
sional analysis  may  be  subdivided  into  three  classes : 

1)  . Results  which  apply  to  problems  requiring  only  C*  continuity 

(C°  displacement  fields,  problems  in  plane  elasticity) 

2) .  Results  which  apply  to  problems  requiring  both  C®  and  conti- 

nuity (coupled  C®  and  displacement  fields) 

3) .  Results  which  apply  to  problems  requiring  only  continuity 

(C^  displacement  fields,  plate  bending  problems) 

2.1  C*  Displacement  Fields  (plane  elasticity) 

It  has  been  shown  in  (6,  7,  8,  9]  that  elemental  arrays  may  be  effi- 
ciently generated  through  the  use  of  ’’precomputed''  arrays  - - that  is  arrays 
which  are  computed  once,  stored  on  permanent  file,  and  then  reused  in  all 
subsequent  applications  of  the  program.  The  new  work  done  by  the  principal 
investigator  and  his  collaborators  has  two  objectives;  the  first  is  to 
show  how  the  hierarchal  C®  elements  (described  in  1.2)  for  a quadratic  functional 
nay  be  formulated  using  precomputed  arrays  thus  yielding  a finite  element 


! 

I 

technique  which  is  especially  suited  to  problems  with  local  rapid  variation 
of  the  function  to  be  approximated.  In  particular,  formulas  for  two-dimen- 
sional (hierarchal)  element  arrays  for  ^u:bitrary  polynomial  order  are  derived, 
based  on  precomputed  arrays.  The  second  objective  is  to  apply  the  combined 
approach  of  hierarchal  elements  and  precomputed  arrays  to  decide  if  a computed 
result  has  "converged".  A common  practice  in  finite  element  analysis  is  to 
solve  a problem  several  times  using  Successively  required  meshes  i.e.  to 
apply  the  procedure  for  h-convergence.  If  successive  analyses  agree  then  it 
is  usually  assumed  that  the  finite  element  approximation  is  accurate.  This 
procedure  can  be  computationally  expensive  when  several  highly  refined  meshes 
are  used.  An  alternative  procedure  is  to  use  p-convergence,  which,  as 
pointed  out  in  1.1,  has  a faster  rate  of  convergence  to  the  true  displacements. 
The  computational  effectiveness  of  the  p-convergence  procedure  is  demonstrated 
numerically  using  hierarchal  elements  and  precomputed  arrays. 

Detailed  formulas  are  given  for  calculation  of  stiffness  matrices  and 
for  calculation  of  polynomial  coefficients  from  nodal  variables.  Hierar- 
chal nodal  variables  are  presented  together  with  some  of  the  favorable 
consequences  of  using  hierarchal  nodal  variables.  Computation  times  for 
stiffness  matrices  are  given  in  terms  of  equivalent  time  units  (e.t.u.)  for 
different  methods.  In  Table  1 average  Central  Processor  Unit  (CPU)  times 
are  given  for  computing  element  stiffness  matrices  for  several  different 
problem  types  and  polynomial  orders  p.  For  comparison  purposes  the  compu- 
tation time  for  an  eight  degree  of  freedom  isopareunetric  quadrilateral  using 
the  general  structural  analysis  program  SAPIV  of  Bathe  et  al  [10]  has  been 
included  in  the  table.  The  time  required  to  generate  this  8x3  matrix  is 
seen  to  be  comparable  to  that  required  for  the  45  x 45  hierarchal  matrix  of 


the  torsion  problem.  This  means  that  the  constraint  method  allows  many  more 
degrees  of  freedom  for  the  same  computer  cost. 


various  numerical  examples  are  analyzed  in  detail.  One  such  example 
IS  a Cantilever  Beam  wit.h  Uniform  Load,  which  we  now  describe. 

2.1.1  An  Example;  Cantilever  Beam  with  Uniform  Load 

It  is  Itnown  that  in  finite  element  problems  where  roundoff  error  is 


significant,  an  Increase  In  the  polynomial  order  may  be  advantageous 
[11].  This  suggests  that  employing  hlerarchal  elements  with  p-convergence 
may  oaks  It  possible  to  use  "poorly-shaped"  elements  - that  Is,  elemants 
of  such  proportions  as  would  lead  to  ill-conditioned  equations  If  low 
order  polynomials  were  used.  To  test  this  Idea,  we  consider  the  plans 
stress  cantilever  beam  shown  In  Fig.  4,  loaded  with  the  edge  stresses 
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The  beam  is  constrained  against  rigid  body  motion  through  the  boundary 

conditions  u - v ■ 0 for  x ■ y ■ 0,  and  u ■ 0 for  x ■ 0,  y ■ -d/2.  The 
five-element  mesh  employed  In  the  analysis  Is  shown  In  Fig.  4.  Table  2 

contal.ns  the  displacements  (In  Inches)  at  point  A and  the  potential 
energy  « calculated  for  the  case  of  a deep  beam  (L  " 2.0  In.,  d " 1.0  In.) 
with  elastic  modulus  of  2.8  x 10^  psl,  and  a slender  beam  (L  ■ l.S  In., 
d - 0.1  In.)  with  elastic  modulus  of  2.8  « 10^  psl.  Pnlssca's  ratio  wae 
taken  as  0.3  and  the  thickness  as  unity.  The  hlerarchal  elements  exhibit 
convergence  In  both  cases,  even  though  the  ratio  of  height  to  width  la 
1:30  for  two  of  the  elements  In  the  mesh  for  the  slender  beam.  Some 
numerical  difficulty  did  occur,  however,  when  an  attempt  was  mada  to 


solve  Che  slender  beam  problem  for  values  of  L ~ 15  In.  and  d ~ 1.0  in. 

Since  a simple  scaling  of  chese  dimensions  led  Co  accurace  resulca,  ic 
may  be  chac  Che  difficulcy  was  noc  due  co  numerical  ill-condiclonlng, 
but  inscead  was  caused  by  Che  Cechnique  used  Co  impose  Che  boundary 
condicions  - chac  is,  Che  arciflclal  sClffness  coefflcienCa  Introduced 
on  Che  diagonal  of  Che  global  stiffness  macclx  were  not  large  enough. 

For  comparison  purposes,  Table  2 also  contains  values  calculated 
by  Mason  and  Case  (12]  using  a fine  mesh  of  low  order  elements.  Since 
no  analytical  solution  exists  for  the  beam  under  the  given  loading, 
these  authors  validated  their  finite  element  model  by  applying  it  to  a 
slightly  different  problem  - involving  Che  same  beam,  but  a different 
loading  - for  which  an  analytical  solution  exists.  The  finite  element 
and  theoretical  results  for  the  displacements  were  found  to  agree  to 
within  0.4S.  In  the  present  example,  the  results  from  [12]  and  from 
the  hlerarchal  element  solution  were  also  found  to  .agree  to  within  0.4Z, 

The  hierarchal  element  results  may  be  considered  closer  to  the  exact 
values,  if  the  good  convergence  of  the  computed  values  observed  in  Table 
2 can  be  taken  as  an  indicator  of  accuracy. 

Table  3 contains  CPU  times  for  the  equation  solver  for  Che  hlerarchal 
element  approach. 

2.2  Couoled  C*  and  displaceiwnt  fields 

In  this  work  the  results  of  Kratochvil  et  al  in  (6)  are  generalized 
to  problems  with  three  independent  displacement  fields.  An  essential  aspect 
of  this  approach  is  to  transform  a triangular  element  T in  the  x-y  plane 
into  a standard  triangle  T with  vertices  at  the  origin  and  at  a unit  distance 
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along  the  horizontal  and  vertical  axes.  Such  a transformation  is  shown 
below  where  C and  n represent  coordinates  in  the  plane  of  the  standard 
triangle  and  x and  y are  local  coordinates  for  the  element  with  the  first 
vertex  of  the  element  coinciding  with  the  origin.  The  other  two  vertices 
and  also  the  three  edges  are  numbered  in  counter-clockwise  order  as  shown. 
The  justification  for  using  the  standard  triangle  is  that  integrations  and 
matrix  inversions  are  performed  with  respect  to  the  standard  triangle.  Thue 


they  need  be  done  only  once  and  the  results  are  stored  and  then  used  in 
all  future  applications  of  the  program.  Computation  of  the  element  stiffness 
matrix  is  thus  reduced  to  computing  a linear  combination  of  a small  number 
of  precomputed  matrices  followed  by  pre-  and  post-multiplication  by  block 
diagonal  matrices.  The  number  of  precomputed  matrices  which  must  be  stored 
are  considerably  reduced  by  choosing  hierarchal  nodal  variables. 
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2.3  C displacement  fields  (place  bending  probl«ims) 

It  was  suggested  by  Dr.  V.  B.  Venkayya  (of  t.he  Analysis  and  Optimi- 
zation Group,  Structures  Division,  Air  Force  Flight  Dynamics  Laboratory, 
Wright-Patterson  Air  Force  Base)  in  a letter  dated  12  December  1975  t.hat 
further  development  of  plate  bending  elements  would  be  useful  in  order  to 
realize  the  full  potential  of  the  constraint  method.  Accordingly,  a 
sophisticated  plate  bending  element,  incorporating  a complete  p^  order 
polynomial  with  p^5  and  corrective  rational  functions  is  now  being 
programmed  and  tested.  We  now  describe  some  of  this  work. 

It  is  well  known  (see  (13],  for  example)  that  exactly  conforming 
(even  at  vertices)  displacement  fields  cannot  be  formed  merely  by  freely 
assembli.ng  finite  ele-oents.  Thera  are  certain  additional  constraint  equa- 
tions which  must  be  satisfied  at  vertices.  The  sLmplest  form  of  these 
constraint  equations  has  been  given  by  Peano  in  (14],  where  a specially 
devised  asse.mbly  procedure  is  presented  which  automatically  enforces  the 
constrai.nt  equations.  Alternative  methods  for  enforcing  exact  conformity 
when  using  displacement  fields  of  arbitrary  polynomial  order  p are: 

a)  . creation  of  super-elements  (or  macro -elements)  of  arbitrary  order 

p^5.  In  these  super-elements  constraints  are  satisfied  withi.n 
each  super-element  leaving  nodal  variables  on  the  boundary  to  be 
freely  assembled  (14]. 

b)  . Use  of  penalty  functions  to  enforce  constraints  [14] 

c) .  Supplementing  pzll  order  polynomials  with  newly  constructed 

corrective  rational  functions.  This  destroys  the  analytic 
character  of  the  approximation  at  the  vertices  but  permits  free 
assembly  of  elements  without  enforcing  constraints  (14,  15). 
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The  last  alternative  is  the  one  which  seems  most  promising.  An  algorit.hm 
has  been  developed  for  a (exactly)  conforming  trianguleir  element 
which  contains  complete  polynomials  of  order  p^5  and  rational  corrective 
functions.  A typical  element,  of  order  p >•  5,  together  with  nodal  variaibles 
is  shown  below. 


For  the  quintic  element  there  are  24  independent  nodal  variables.  The 
shape  functions  for  each  nodal  variables  have  been  given  explicitly  in 
[14,  15]  in  terms  of  triangular  coordinates.  The  shape  functions  corres- 
ponding to  second  order  tangential-normal  derivatives  at  the  vertices  are 
the  rational  functions  given  below;  all  other  shape  functions  are  poly- 


nomials. 
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It:  ia  ia)?ortanc  to  observe  that  although  rational  functions  are  usee 
in  Che  basis,  all  terms  which  appear  in  Che  elemental  stiffness  matrix 
can  be  integrated  explicitly  without  recourse  to  numerical  quadrature. 
This  was  proved  in  [16).  An  algorithjn  based  on  a hierarchical  family  of 
elements  using  corrective  rational  functions  has  been  programmed  and 
IS  now  being  tested  on  numerical  examoles. 
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In  addition  work  ia  proceeding  on  the  development  of  three  dimensional 
finite  elements  for  use  in  the  constraint  method  as  follows. 

2 . 4 .gpdal  Variables  and  Shape  functions  for  a hierarchical  family 
of  solid  tetrahedronal  elements 

A table  of  basis  functions  for  a triangular  (two  dimensiona  ) 
hierarchical  family  was  given  in  [li] . This  is  now  being  generalized  to 
a table  of  basis  functions  for  a tetrahedronal  (three  dimensional) 
hierarchical  family.  From  this  new  table  it  will  be  possible  to  generate 
shape  functions  for  solid  C*  elements  and  (if  needed)  corresponding  nodal 
variables . 

2.5  Parers  craoared  for  oublication  and  cresentation  at  conferences 


The  work  described  in  2. 1-2. 4 will  be  reported  to  the  professional 
community  in  the  following  papers  to  be  submitted  for  publication  in 
journals . 

1) .  "Hitrarchal  Finite  Elements  and  Precomputed  Arrays",  by 

Mark  P.  Rossow  a.nd  I.  Norman  Katz,  (to  be  submitted  to 
Int.  J.  for  Num.  Method  in  Engr.). 

2)  . "Modal  Variables  for  Conforming  Finite  Elements  of  Arbitrary 

Polynomial  Order",  by  I.  Norman  Katz  and  .Mark  P.  Rossow,  (to 
be  submitted  to  Computers  and  Mathe-matics,  with  Applications). 

3) .  "A  family  of  C^  triangular  elements  containing  complete  poly- 

nomials of  arbitrary  order,  for  application  to  problems  in  plate 
bending",  I.  Norman  Katz,  Darna  A.  Szabo  and  Olive  Liu,  (in 
preparation) . 

The  following  papers  have  been  accepted  for  presentation  at  conferences. 

4) .  "Hierarchical  Approximation  in  Finite  Element  Analysis",  by 

I.  NOrman  Katz,  International  Symposium  on  Innovative  Nu.merical 
Analysis  in  Applied  Engineering  Science,  Versailles,  France, 

May  23  - 27,  1977. 

5) .  "Efficient  Generation  of  Hierarchal  Finite  Elements  Through  the 

Use  of  Precomputed  Arrays",  by  M.  P.  Rossow  and  I.  N.  Katz, 
Second  Annual  ASCE  Engineering  .Mechanics  Division  Specialty 
Conference,  North  Carolina  State  University,  Raleigh,  NC 
May  23  - 25,  1977. 
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6).  "C^  Triangular  Elements  of  Arbitrary  Polynomial  Order  Containing 

Corrective  Rational  Functions",  by  I.  Norman  Katz,  SIAM  1977 
National  Meeting,  Philadelphia,  PA,  June  13  -15,  1977, 


Table  1 


Computation  Tine  for  Element  Stiffness  Matrices 


Problem 

Type 

Plelds 

P 

Size  of 
Matrices 

CPU  Time 
In  e.t.u 

Plane  stress 
(energy  method) 

Two  displacements 

5 

42 

X 42 

2.1 

Torsion 

(energy  method) 

Prandtl  stress 
function 

8 

45 

X 45 

1.3 

Plane  stress 
(least  squares) 

Two  displacements 
and  three  stresses 

3 

50 

x-50 

3.4 

Plane  stress 
(energy  method) 
SAPIV  (10] 

Two  displacements 

1“ 

8 

X 8 

1.2 

Plate  bending 
(energy  method) 
Ref.  7 

One  displacement 

5^ 

18 

X 18 

0.9*^ 

!l 

i 1 


a Bilinear 

b Incomplete  fifth  order  polynomial 

c Matrix  multiplication  routine  used  In  defining  an  e.t.u.  was  not 
specified 
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TabU  : 

Raauica  tor  Cantilever  Bean  Problem 


P 

DOF 

v^.IqS 

3 

59 

-8.683437 

4 

99 

-8.757509 

5 

149 

-8.760231 

6 

209 

-8.764954 

7 

279 

-8.764951 

8 

359 

-8.765809 

1152 

Bilinear 

Ref. 

elemenca 

12 

2447 

-8.7494 

2160 

Bilinear 

Ref. 

el ements 
12 

4703 

- 

, d-1  L-1.5, 


Tt.lO® 

w.lO^ 

-3.6164721 

-4.027224 

-1.7374483 

-3.6227009 

-4.033644 

-1.8010797 

-3.6227475 

-4.033650 

-1.8010797 

-3.6227711 

-4.033656 

-1.8010798 

-3.6227734 

-3.6227744 

-4.0187 

m 
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Tabla  i 

CPU  Tla*  iUqulrad  to  Sglv  Equaclona 
for  Oaap  Cantllavar  Baaa  Problem 

CPU  Tina 


p 

"l 

la  a. 

3 

62 

20 

4 

102 

36 

5 

132 

36 

6 

212 

96 

7 

282 

174 

8 

362 

298 
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On  September  3,  1977,  the  Principal  Investigator,  I.  Norman 
Katz,  Jointly  with  HarnA  A.  Szabo  presented  a seminar  on  their 
current  results  at  the  Air  Force  Flight  Dynamics  Laboratory,  Wrlght- 
Patterson  Air  Force  Base.  The  title  ot  the  talk  was: 

"Advanced  Stress  Analysis  Technology" 

An  abstract  of  the  talk  is  enclosed  with  this  report.  Forty  to 
fifty  people  attended  the  seminar. 

The  prlncip<il  contact  at  Wrlght-Patterson  was  Dr.  V.  N.  Ven- 
kayya.  Aerospace  Engineer,  Analysis  and  Optimization  Branch,  Struc- 
tures Division.  Other  people  with  whom  we  spoke  personally  were:  Dr. 
N.  3.  Khot,  Dr.  L.  Berke,  Dr.  D.  W.  Quinn  and  Mr.  N.  D.  Wolf  of  the 
Analysis  and  Optimization  Branch;  Dr.  J.  Oallagher  and  Mr.  R.  M.  Bader 
of  the  Structural  Integrity  Branch. 
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A3STRACT 

With  one  exception,  all  finite  elexent  soft^.'are  systens  have  eleneat 
libraries  in  which  the  approxination  properties  of  eleaents  are  frocen. 

The  user  controls  only  the  nunber  and  distribution  of  finite  elesents. 

The  exception  is  an  experlnental  sofrware  systen,  developed  at  Washington 
University.  This  systea,  called  COMET-X,  employs  conforming  elements 
based  on  complete  polynomials  of  arbitrary  order.  The  elements  are  hierar- 
chic, t.e.  the  stiffness  matrix  of  each  element  is  embedded  in  the  stiff- 
ness matrices  of  all  higher  order  elements  of  the  same  kind.  The  user  con- 
trols not  only  the  number  and  distribution  of  finite  elements  but  their 
approximation  properties  as  well.  Thus  convergence  can  be  achieved  on 
ficted  mesh.  This  provides  for  very  efficient  and  highly  accurate  approxi- 
mation and  a nev  method  for  computing  stress  intensity  factors  in  linear 
elastic  fracture  mechanics.  The  theoretical  developments  are  outlined, 
numerical  e.xamples  are  given  and  the  concept  of  an  advanced  self-adaptive 
finite  element  software  systea  is  presented. 


